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SECTION 1 


INTRODUCTION 
1.1 BACKGROUND 


The Experimental Prototype Automatic Meteorological System (EPAMS) is 
part of an overall development program of the U, S. Army to provide a centralized 
data system for use in supporting future operations of the Field Army. Under the 
EPAMS development program, the Atmospheric Sciences Laboratory, White Sands 
Missile Range (WSMR) required implementation of an automated procedure for esti-~ 
mating the cloud trajectory and dispersion characteristics downwind from sources 
of contaminate material released at any arbitrary point within the atmospheric mix-~ 
ing layer. WSMR wished to implement these automated procedures using available 


technical approaches to prevent undue delays in the EPAMS development schedule. 


The technical problem of estimating wind fields over complex terrain has 
been under continued investigation at WSMR. Because of the computer space/time 
criteria established for EPAMS, the model selected for use in EPAMS had to be 
relatively simple in concept. The computer algorithms developed by Tingle and 
Bjorklund (1973) for WSMR's use in predicting wind fields in complex terrain were 
ideally suited for this purpose. The algorithms, based on the shallow-fluid equa- 
tions of oceanography, have simple input requirements, reach stable solutions in a 
relatively few iterations and thus conserve computation time, and have been shown 
to produce wind field estimates agreeing with observed estimates in complex ter- 
rain surrounding WSMR and Dugway Proving Ground, Utah. Similar criteria applied 
to the selection of dispersion models for use in EPAMS. The generalized dispersion 
models developed under previous Government contracts for predicting concentration, 
dosage and deposition downwind from nearly-instantaneous and continuous volume 
sources (Cramer, et al., 1972; Dumbauld, et al., 1970, 1973) had been validated 


for use in conjunction with military activities. These computerized models were 
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developed with the intention of providing the U. S. Army with a universally adaptable 
means of calculating dispersion characteristics for all types of terrain, meteorolog- 


ical conditions and types of sources. 


Thus, the technical data base requisite for development of an automated pro- 
cedure for predicting cloud trajectories and dispersion characteristics in complex 
terrain existed, The principle technical objectives of the work described below were 


to utilize the existing technical data base in the: 


e Development of a diagnostic/prognostic analysis routine 
for evaluating the depth and wind velocity of the atmo- 


spheric mixing layer in association with complex terrain 


e Development of a diagnostic/prognostic application rou- 
tine for evaluating the transport and diffusion of contam- 
inant material released from a source within the mixing 


layer 


and incorporation of these routines into EPAMS as subordinate elements. 


Figure 1-1 is a simplified schematic diagram of the EPAMS functional ele- 
ments illustrating the interface of the Analysis and Application Routines or Phases 
with EPAMS. The Applications Phase is called by the EPAMS Executive as a result 
of a request for information routed through a Message Management Routine. The 
Applications Phase extracts the necessary meteorological and source information 
from the Data Base Management Routine and calculates the transport and diffusion of 
contaminant material released to the atmosphere. The solution is output to the Data 
Base, where it is made available to the EPAMS Executive for forming a response to 
the information request. The Analysis Phase is called by the EPAMS Executive as 


a result of a standing request for a routine analysis of meteorological conditions in the 
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surface mixing layer. The Analysis Phase extracts upper-air and surface meteorolo- 
gical data, including U, S. A. F. Global Weather Central (GWC) "fine mesh" pre- 
diction data, from the Data Base. As shown in Figure 1~1, the Analysis Phase is 
comprised of two major subelements, a mixing-layer analysis routine and mesoscale 
wind-field model. The mixing-layer analysis routine calculates initialization param- 
eters for use in the mesoscale wind-field model as well as additional meteorological 
inputs for direct use in the Applications Routine. The mesoscale model calculates 
wind-field patterns and mixing-depth heights for use in the trajectory-transport sub- 
element of the Applications Phase of EPAMS. The Analysis Phase does not commun- 
icate directly with the Applications Phase, but supplies solutions to the data base for 


subsequent use as required by the Applications Phase. 
1.2 ORGANIZATION OF THE REPORT 


A description of the mesoscale wind-field model and the transport diffusion 
models respectively used in the Analysis and Application Phases is provided in Sec- 
tion 2 below. Section 3 contains a description of procedures used in the Analysis 
Phase to develop meteorological inputs to the mesoscale wind-field model and for 
direct input to the diffusion models from meteorological data extracted from the 
EPAMS data base. Details of the computer program structure of theAnalysis and 

Application Phases are given in Section 4, Section 5 contains example cases illus- 
trating the use of the computer program. Detailed user instructions for the com- 
puter program routines are given in Appendix A and Appendix B contains a listing 


of the computer program routines. 
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SECTION 2 


MESOSCALE WIND-FIELD AND 
TRANSPORT/DIF FUSION 
MODELS 


2.1 MESOSCALE WIND-PIELD MODEL 


The Analysis Phase developed for EPAMS relies on 4 numerical modeling 
technique developed for WSMR by Tingle and Bjorklund (1973) to calculate details 
of the mesoscale wind-field and depth of the surface-mixing layer over complex 
terrain, The computer algorithm, based on the shallow-f{luid equations of ocean- 


ography, describes the motion of fluid in two dimensions by the expressions: 


du ou ou o) 
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ot ox oy is i) 
where 
u,v = x, y components of the fluid velocity in the layer next to 
the terrain 
Hn = depth of the surface mixing-layer 
Zn = height of the terrain 
g' = reduced acceleration of gravity 
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g = acceleration of gravity 
Po = density of fluid in the surface mixing layer 
mo = density of fluid in the deep passive layer above the 


fluid in the surface mixing-layer 


Thus the equations can be used to describe mesoscale wind-patterns over complex 
terrain when the atmosphere above the terrain can be distinctively divided into at 
least two layers of air of different densities, Such an atmospheric condition arises 
when the mixing layer above the terrain is capped by an elevated inversion. The 
warmer upper layer interacts with the cooler lower layer only in that its lesser 
density reduces the speed of gravity waves introduced into the lower layer. In the 
algorithm, Tingle and Bjorklund have written Equations (2-1), (2-2) and (2-3) in their 
momentum form. Mesoscale wind patterns are generated by impulsively accelerating 
the momentum in the lower layer to a constant value and, after a empirically deter- 
mined model time (see Section 4), the flow in the region of interest is assumed to 


be stabilized (steady state). 


The computer algorithm utilizes an explicit method of the Lax-Wendroff 
type of finite differencing technique using nine grid points and one time level in 
combination with a nine-point low-pass filter. The differencing scheme is designed 
such that the finite differencing forms are centered to preserve high-order accuracy, 
The grid arrangement for this centering scheme and finite differencing form of the 
equations are fully described in the report by Tingle and Bjorklund (1973), As noted 
above, high-frequency oscillations generated in critical flow situations near regions 
where hydraulic jumps occur are removed by a filtering technique. The filtering 
technique used in the algorithm is based on a suggestion from Shapiro (1970) and 


has the advantage that its smoothing effect is easily determined from its response 
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function, In the Analysis Phase, the filter is applied every seventh time step. It is 
well known that the length of the time step At used in obtaining the solution also 
effects the stability of the finite difference solution. ‘The time step is altered 


automatically in this algorithm in accordance with the expression 


ii] 


x ‘ 
Ag ZS eines min 


max 


where the stability parumeter A is dependent on the grid spacing Ax and varies 
between zero and one. For the present application using a minimum spacing (Ax ay, 
of 5 kilometers, the stability parameter has been set equal to 0.4. Because the 

entire solution grid must be searched to determine the maximum wave speed in the 


denominator of Equation (2-4), the value of At is calculated every second time 


step rather than every time step to save computation time. 


A sensitivity analysis of the two-layer mesoscale wind field model cur- 
rently in progress under separate contract with U. S. Army Dugway Proving Ground 
has indicated model instability can frequently occur in localized terrain regions 
where the inversion height Hn is initialized, or becomes adjusted by the model, to 
altitudes within a few tens of meters of the higher terrain elevations. The momen- 
tum conservation properties of the computational scheme are such that wind veloc- 
ities in these localized regions may become unreasonably large as the layer thick- 
ness approaches a minimum value. In order to prevent this type of solution insta- 
bility, the computer algorithm has been redesigned so that the u and v compon- 
ents of the solution at any given time step are automatically set to zero and the 
thickness raised to a minimum depth of 30 meters if the calculated depth for a grid 


point is less than 30 meters. Limited experience in the application of this feature 
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of the mesoscale model indicates that computational stability is preserved without 
undue loss of momentum over the entire grid and with the result that wind fields in 
the immediate vicinity of higher elevations under these conditions appear to be rea- 


sonable (see Example in Section 5), 


The computerized algorithm uses an expanding grid in the boundary area 
beyond the terrain region of interest to damp outgoing waves and eliminate incom- 
ing waves based on the technique suggested by Lavoie (1972), Grid points in the 
boundary region are spaced at 10, 26, 10, 80 and 160 kilometers from all four sides 
of the interior grid. The terrain height at all points in the boundary region are set 
equal to the minimum elevation within the interior grid, Since the shallow-water 
equations are written in their momentum form, the initial momentum is set equal 


to the momentum used in initializing parameters in the interior grid. 


Readers interested in the form of the finite-differenced equations and other 
details on which the mesoscale wind-field computer algorithm is based are referred 


to the paper by Tingle and Bjorklund (1973). 
2.2 TRANSPORT-DIF FUSION MODELS 


beara k Trajecto ry/ Transpo rt Routine 


The Trajectory-Transport Routine incorporated in the Applica- 
tions Phase of EPAMS shown in Figure 1-1 is designed to calculate cloud trajectories 
and mean-layer winds and mixing-layer depths at fixed intervals along the trajector- 
ies using the output from the mesoscale wind-field model described above. The 
trajectory of a cloud released at any point within the model solution grid is calcu- 
lated using a simple Euler predictor-corrector scheme and nine-point interpolation 
procedure. The computer program uses wind vector data from the mesoscale wind- 


field model solution at nine grid points surrounding the source (point of release) and, 
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through the interpolation procedure, estimates the u and v components of the 

wind at the point of release. The first approximation to the first point on the tra- 
jectory curve is then determined by resolving the components and placing a point 

at a fixed interval along the resultant vector, ‘The fixed distance interval is set at 
one-tenth the grid spacing used in the mesoscale wind-field solution matrix so that 

in the present program a fixed interval of 500 meters is used, Wind velocity com- 
ponents at this point are again estimated from the nine-point interpolation proced- 
ure. These new u and v components are then averaged with the previously deter- 
mined components at the starting point. The final estimate of the first point on the 
trajectory is then calculated by resolving the averaged components and moving a 
fixed interval along the resultant vector. This final estimate of the first point on 

the trajectory is then treated in the same manner is the initial starting point (source) 
and the interpolation and iteration scheme repeated to find the second point on the 
trajectory. This procedure is repeated until the trajectory passes outside the grid 
area containing wind-vector data or the number of iterations exceeds a predetermined 


limit. 


The interpolation method is based on an evaluation of a truncated two-dim- 


ensional Taylor's expansion formula (Carnahan, Luther and Wilkes, 1969, p, 430 
y ’ 1 


given by 
he ee 
my) ~ E tht the +e + hke FS a- 
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i = x+=< 
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and fj is either the u or v component at Xo y; , the nearest grid point to the 


point of interpolation (x,y). The partial derivatives f,f,f , f and f are 
S Yo Es xy yy 


estimated by divided differences. 


The depth of the mixing layer at fixed intervals along the trajectory is also 
interpolated from data provided by the mesoscale wind-field model using a four-point 


bivariate interpolation formula given by 


ix ov = on, = + = 
He (x, y) (=p) (CE = ic) ae ee (i=5q) Hh . 
i,j itl, j 


+ q (1 - p) a " r py Be 
i, jr Tp Ba pug 


a: SM ht A A el ak mn a lL ae Sa 


where 
po 8/4) 
7 = Gey Weer y;) 
and Hy F Hs ; t » and Hy are the mixing layer depths at the 
i,j i+1, j i, jt1 iP ajtt 


four surrounding grid points. 


2268 Diffusion Models Routine 


The Diffusion Models Routine shown in Figure 1-1 as an element 
of the Applications Phase of EPAMS is comprised of volume source diffusion models 
for instantaneous and continuous sources and models for use in predicting the buoy- 
ant rise of plumes from continuous sources. These models are used in the Diffusion 


Models Routine to estimate isopleths of dosage and concentration, with or without 
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depletion due to simple decay, gravitational settling and precipitation scavenging, 
* along the cloud trajectory determined by application of the Trajectory-Transport 


Routine. 


Instantaneous Volume-Source Model 


The volume-source model for use in EPAMS is similar to the models for 


instantaneous sources described by Cramer, et al, (1972) and by Dumbauld, | 


Bjorklund and Bowers (1973). The concentration for an instantaneous volume-source | 
can be expressed as the product of four terms: | 
xy = {Peak Term} {Lateral Term} {Vertical Term} (Depletion Term} (2-6) 


The Peak Term describes the concentration at the cloud centerline and is given by 


the expression 


Q 


Peak Term (2-7) 


aes GO 
xn yn zn 


where 


Q = total source strength 
‘. = standard deviation of the concentration distribution in the 
; alongwind direction at the nt) distance along the cloud tra- 
jectory 
0 2 standard deviation of the concentration distribution in the 
4 y crosswind direction at the nth distance along the cloud tra- 
j jectory 


= standard deviation of the concentration distribution in the 
vertical at the n' distance along the cloud trajectory 
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The Vertical Term describes the vertical distribution of material at the 
h ‘ : , 
n distance along the cloud trajectory as modified by reflections of material at 
the top of the surface mixing layer and at the earth's surface and by gravitational 


settling of particulates 


( N Ax 

0 QaH o-H' +2+V_ — 

Vertical _ PS re nhs eile ee nmi “n 
Term 2 (7 


> AX, 2) 
aH Pag sie fas 


at 1 cpa mee 
{ y 2 LN a 5 
ie Z,n 
-8) 
( N 
22a HL tH' -2-V 
ae ele 
+ > <r expl- — Oe TES eee Corer son” | -a ES.. 
2 u 
1 ] Zn 
2aH ed SS Sale 
ps a aia hacer aerate ee 
P 2 o 
ZN 
F is ee 0 
| where, for convenience in writing the Vertical Term, 0 (zero to the zeroeth power) 


‘ is set equal to unity and 


| 
| 
i 
' 


He = effective release height (see Equation (2-20) ) 
z = ealeulation height 
th 
Nn - = mixing layer depth from mesoscale model for n 
’ 


distance along the trajectory 


r = fraction of material reflected at the surface (1 for 
complete reflection and 0 for no reflection) 
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vi = settling velocity for material in a given size cate- 
gory 
Ax, = incremental distance along the cloud trajectory 
a = mean wind specd in the layer containing the cloud 


at the nth point along the trajectory (see Equation 
(2-17) below) 


The Lateral Term is given by 


1 
Lateral Term = exp wr os (2-9) 
i yn 
where 
y =  crosswind distance from the centerline of the cloud 


trajectory 


The Depletion Term refers to the loss of material by simple decay processes 
and by precipitation scavenging. The form of the Depletion Term for each of these 


processes is: 


| Y AX 
, (Decay) exp |-k > — (2-10) 
m1 “n 
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(Precipitation - ae 


Scavenging) en ie Git 


whe re 


k = decay coefficient or fraction of material lost per 


unit time 


A washout coefficient or fraction of material removed 


by scavenging per unit time 


The subset of equations describing the standard deviations of the alongwind 


go, crosswind 
xn 


xn 


where 
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and vertical o concentration distributions is: 
VAN 


Ln tes 0.6 Au 
Hy) + y = ie AX y 
( Ss = ? : 
1.3 x(n=-1) 1.3 u n x(n-1) 


alonewind dimension of the cloud 


wind speed shear at the nth point along the tra- 
jectory 


u 

5,n (abn ait ) 
ane Tern een 

: 3 sete th 
wind speed at a height of five meters at the n 
point along the trajectory 


Fee Ol ce 
wind power-law coefficient at the n point along 
the trajectory 


top of the cloud at the nth point 


\ H' +2,15¢ 3; 2, < H 
am ron m,n 
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base of the cloud at the nth point 
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(2-14) 
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= Pp 


a (B. . =o) 8 * rw.) 


a Pn) A ae Pn) 


m,n 


mean layer wind speed from mesoscale 
model at nth point 


mean wind speed in the layer containing 
the cloud 


5, Tam Bn , 


p 
(1 +p) @) 9 ( 


a, (x Pn) - 20 * Pn) 


lin Hina) 


incremental distance between the n points 
equal to 500 meters 


standard deviation of the wind aximuth angle 


at the nth point 
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(2-17) 


(2-18) 


(2-19) 


(2-20) 
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T = source function time ~ 2.5 seconds for instan- 
taneous sources 


on = ten-minute standard deviation of the wind azi- 
muth angle 


6 \uayP 

"on i rR n ms oe = B23} 

’ | . E, oO 

where 
n 

u = 7/180 AD 2-22 
oe On /180) (H' /5) ( ) 
oF ‘ = standard deviation of the wind elevation angle 

’ 


The cloud expansion coefficients a (lateral) and 6 (vertical) in the above equations 


are set to unity for instantaneous sources, 


The concentration Xn is calculated under the additional restriction that 


< » JE > the values of o and o are increased proportion- 

Xn = yet Xa 7 Xe? yn x,n sa 
ately so that Ni = Xn 1 and the calculations continued. 

Dosage for an instantaneous source is obtained from the expression 

Xn ¥ ee 
Dye ee = : (2-23) 
n u 
Zn 

where 

u = mean wind speed at the calculation height z and 


nth point on the trajectory 
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Continuous Source Model 


The continuous-source model for use in EPAMS is a modified version of 
the Gaussian model for continuous sources described by Pasquill (1962) in which 


the concentration is also expressed as the product of four terms: 


Xo - = {Peak Term} {Lateral Term} {Vertical Term} {Depletion Term} (2-25) 
, 


The Peak Term of the continuous source concentration model is given by the expres- 


sion 
Q' 
Peak Tenn: = sae (2-26) 
Ir Oo. Oo Ul 

zg yn on 
1 
] 
} where 

Q' = source strength expressed as a rate 

ot = standard deviation of the concentration distribu- 


tion in the crosswind direction for a continuous 
source at the nth distance along the plume tra- 
jectory 


AIS ee lt A 


and uw and 7, ate respectively defined by Equations (2-17) and (2-19), The Ver- 


tical, Depletion and Lateral Terms (with On replaced by Si where appropriate) 
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are identical to the similar terms described for the instantaneous volume source. 


The subset of equations describing o' is: 
yr 
y oO 
te fa 
(0) 
a! at x “aed (2-27) 
yn Ac,n n fee | 
pp. 
where 4 
\ m 
o te 
Ac,n A,n\180/ \5 } (2-28) 
and the remaining parameters have previously defined. The default values of a 
and 8 used in the program for continuous sources are @ =0.9 and B= 1, 
Plume-Rise Models 
The effective height H of a buoyant plume is given by the sum of the phys- 
ical stack height h and the buoyant rise Ah. For a neutral or unstable atmosphere, 
the effective stack height is given by (after Briggs, 1971; 1972) 
1/3 
\is.os | 2 T.\y Sis 
He =) hee i= ow ee h 4 (2-29) 
| uihys s I \ 
eC 
whe re 
- = the inner radius of the stack at the exit (m) 
18 
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3 -1 
w = the stack exit velocity Paine ) 


¥ = the stack exit temperature (°K) for the continuous 
(s 
source 
T. = ambient air temperature (°K) ; 
Ss 
f = empirical correction factor 
h = actual stack height 
ufh} = mean wind speed at stack height 
0) 
ses 
so (2) 2-30 
2, NALS) ( ) 


The factor f, which limits the plume rise as ath} approaches the stack exit veloc- 


ity w, is defined by 


| 1 ; uth} < w/1.5 


3w - 30 fh} t 
ee) ; w/1.5 < uth} <w (2-31) 


0 ; ath} = W | 


Equation (2-29) is used in the EPAMS program when the net radiation index 
equals or exceeds 2, U {h} equals or exceeds 5 meters per second, or the potential 
temperature lapse rate b is equal to or less than 0. When these conditions are not 
met, the program uses the corresponding Briggs (1971) rise formula for a stable 


atmosphere given by 


EG 


RSE ER SRR 


h + 2.397 {| =-r 
uihy < 


(wr AGE ; 
f = j= COs |(=—=7- 
libs & 


u ihy 


where 

stability parameter 
6 2-33 
7 (2-¢ ) 
s 


height-weighted mean potential temperature lapse 
rate over an approximately 200-meter interval 
above h 


z~ht200 
> By 
Bu 
zh 


Z~h+t200 


potential temperature at height Z 
Zz 


potential temperature at height Z +] 


Selb SIE RR AN ERE RB I RN 0 RR 


TSS aR TART 


SECTION 3 
DATA ANALYSIS PROCEDURES 


The EPAMS data base contains data blocks of rawinsonde, USAF GWC pre~ 
dicted wind and temperature data, and surface meteorological information, Section 
3 describes the automated data analysis procedures required to utilize this data in 
the development of parametric inputs to the mesoscale wind-field and diffusion models. 
The Mixing-Layer Analysis Routine described in Section 3.1 below is designed to pro- 


vide initialization parameters required by the mesoscale wind-field model. The Dif- 
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fusion Model Routine, described in Section 3.2, uses surface meteorological data and 
mixing layer parameters provided by the mesoscale wind-field model as direct inputs 


to the diffusion model and for use in defining other model input parameters. 


3.1 MIXING-LAYER ANALYSIS ROUTINE 


a EN eh ffm om 


: 
} 
: 


The Mixing-Layer Analysis Routine uses a relatively simplistic approach 
to utilize the various types of data available for the EPAMS data base in the calcula- 
tion of initialization parameters for the mesoscale wind-field model. The decision 
to maintain a simple approach in this initial effort to develop a comprehensive pre- 
diction program was made so that the effects of logic decisions on the derived param- 
eters could be easily traced and analyzed. Based on this reasoning, the analysis 
routine is currently designed to perform a modified "single-station analysis” rather 
than a more complex analysis of the pressure, wind and temperature fields over 
the entire region of interest. It should also be noted that some of the logic decisions 


made in accepting, rejecting and analyzing the data are somewhat arbitrary at this 


maaan 


a 


point in the development of the program. Where time has permitted, we have 


attempted to study the effects of the various decision processes. However, until 
the program is used to predict dosage and concentration fields under a wide variety 
of meteorological conditions, the total implications of many of the decision processes 


will remain unknown. 


ie 


| 
| 
{ 
| 


Features of the Mixing-Layer Analysis Routine are portrayed in the 
schematic logic diagram shown in Figure 3-1. Inspection of Figure 3-1 shows 


that the key features of the procedure are: 


Rawinsonde data from the closest station with a radius 
of less than 100 kilometers from the release point and 
less than three-hours old at the time of request will be 
used without modification to obtain estimates of the 


input parameters to the mesoscale model 


Rawinsonde data measured at the closest station less 
than 100 kilometers from the release point, greater 
than three-hours old and less than twelve-hours old 

are subject to modification using GWC predicted and 


surface winds and temperatures 


Climatological estimates of mixing-layer depths and 
surface wind speeds are used as mesoscale model 
inputs if all rawinsonde data is more than twelve hours 
old or if there are no stations within a 100 kilometer 


radius of the release point 


Climatological estimates are also used if the rawinsonde 
data is greater than six-hours old and the rawinsonde 
surface temperature and 850 millibar temperatures 

have not been modified because of unavailable or 


inadequate GWC and surface data 


GWC data is used to modify the rawinsonde temperatures 


and winds at the 700 and 850 millibar levels; rawinsonde 
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Figure 3-1, Schematic diagram illustrating procedures for obtaining inputs to the 
mesoscale wind field model. 
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temperatures and winds at the 700 millibar are replaced 
if the GWC predicted temperature for the hour of 
interest is within 5 degrees Celsuis; rawinsonde tem- 
peratures and winds at 850 millibars are replaced 

if the temperature from the reporting surface station 
nearest the time and point of release is within 5 
degrees Celsuis of the GWC predicted temperature 


at the 850 millibar level for the hour of interest 


e Surface temperature and wind data from the closest 
station less than 50 kilometers from the release 
point, less than three-hours old, and with a sur- 
face pressure within 10 millibars of the surface 
pressure from the rawinsonde observation, is used 
to modify the surface rawinsonde observation when 


it is more than three-hours old 


e Surface wind speeds from the closest station to the 
point of release are used as default values for the 
mean-layer wind speed if the rawinsonde data is 


not analyzed 


Climatological estimates of the mixing layer depth shown in Table 3-1 are used as 
default values where the rawinsonde data is not available or is not analyzed for 
some reason, These estimates of the median mixing depths as a function of month 
and time~-of-day are based on the work of Norton and Hoidale (1975), who analyzed 

] 8, 236 rawinsonde observations made at WSMR between 1961 and 1972 using the tech- 


nique described by Holzworth (1967). 


As shown in Figure 3-1, the rawinsonde is analyzed to determine the 


mixing-layer depth, mean-layer wind direction and speed and reduced gravity 
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factor. In developing an automated procedure for determining the mixing-layer 
depth, we studied the use of four methods for establishing that a surface-based or 
elevated inversion was indicated by the rawinsonde observations. The simple cri- 


teria associated with the four methods for determining the base of an inversion 


were: 
Method Criteria 
‘ Ae > 01 
Az 
2 aS ep 
VANS ms 
Ad. 
3 Sa eM) 
EZ ae 
AT 
4 SS =e 10° 
ZONY 5 
where 
Ad A 
yen lapse rate of potential temperature 
AL lapse rate of temperature 
{ Az 
Ad. 
a lapse rate of wet-bulb potential temperature 
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= lapse rate of virtual temperature 


The computer program incorporating these criteria was used to analyze twice daily 
rawinsonde observations from Salt Lake City, Utah for the first five days of each 
month during the period June 1972 through May 1973 and for all the month of January 
1973, The heights of the bases of stable layers obtained from the computer program 
were then qualitatively compared with those selected by meteorologists from plots of 
the data. All the criteria produced bases of stable layers in agreement with those 
selected by the meteorologists for cases in which clearly-defined stable layers 

were present, In our view, however, the criteria based on the lapse rate of virtual 
temperature (Method 4) performed more consistently than the other critcria when 
stable layers were not as clearly defined and for all seasons of the year. For this 


reason, the procedure finally used in the program is hased on this criteria. 


The program accepts the raw or modified rawinsonde data and begins the 


selection process by checking 


BE ce 
v,n = v,n V; n (3-1) 
Az Z ue 
n nt+1 n 


in the first height interval above the surface. If the criterion for the base of a stable 
layer is met e.g. Pei ae iA Z, 2-5 x 10-4) and the height interval Az is 
greater than or equal to 100 meters, the program defines a surface-based inversion 
and sets, by default, the surface mixing depth to 30 meters. If Az, is less than 


S= 10-4 ‘heck / 
100 meters but ev, /A 2,2 5 x 107%, the program then checks si 8 9/825: 


>- =< € + > P 
If AT, 5/42 >-5x10™ and (Az, Az,) > 100 meters, the program checks 
the quantity 


: 


} 


: Se! T ~ 
(at) - My nte Van po) Nae al (3-2) 
ae ‘nian a 


If (AT/Az) > —5 x to", the program again defines a surface-based inversion. 
If no surface based inversion is found, the program increments n and proceeds 
through the same operations until the base of an elevated inversion is found or until 
Zed exceeds 3000 meters AGL. If zd exceeds 3000 meters AGL, the surface 
mixing depth is defaulted to 3000 meters. 


After the depth of the surface mixing-layer Hn is established, the program 
analyzes the rawinsonde data between the surface and a to obtain estimates of the 


mean-layer wind components and the reduced gravity factor using the expressions: 


Hy 
zu a (+4 2 - 
“En = = Gets (3-3) 
bere} Gra ., es) 
H 
yy “nh a+ of 
= S17 eine a ae 3-4) 
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where T is the maximum temperature in the stable layer above H . In 
I, max : m 


the program, g!' is restricted so that if it exceeds 0.3 it is set to 0.3 and if less 


than 0.1 is set to 0.1. 
Se2 DIFFUSION MODEL ROUTINE 


The meteorological and source parameters required by the diffusion models 
are given in Table 3-2. Of the meteorological parameters shown in the table, the 
user need only supply the program with values of r and A if they differ from the 
default values of 1 and 0, respectively. The program user must supply all the 
source input parameters to the program and, in addition, denote whether the source 
is instantaneous or continuous. For program test purposes, the program now con- 
tains two default continuous sources corresponding to the American Smelting and 
Refining Company stacks associated with their lead and copper refineries in El 
Paso, Texas. The source inputs for these stacks were supplied by WSMR. Default 
values of ¥.. and k equal to zero are used in the program. 

Since the turbulence parameters ¢! and o' , their power-law coef- 

A,n E,n 
ficients for height dependency m and n, and the power-law coefficient p, are 
not routinely measured at meteorological stations, a means of assigning appropriate 
values to these meteorological parameters was required. Swanson and Cramer 
(1965) made a comprehensive study of the height dependence of go, asa function 
of wind speed and time-of-day using measurements from a 62-meter tower at WSMR, 
The decision was made to use the results of this study in conjunction with a stability 
classification system, similar to the Pasquill definition of stability categories 
(Turner, 1964) used by the Environmental Protection Agency, to select and assign 
values to these parameters. Figure 3~2 is a schematic diagram showing the logic 
sequence used by the program in defining the parameters and their use in the dif- 


fusion models. 
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TABLE 3-2 


METEOROLOGICAL AND SOURCE INPUTS 
FOR THE DIFFUSION MODELS 


(a) Meteorological Inputs 
Net radiation index 


Depth of the surface mixing layer at the nt) point in the 
trajectory 


Mean layer wind speed in the surface mixing layer at the 
nth point on the trajectory 


Standard deviation of the wind azimuth angle at the nth 
point on the trajectory 


Power-law coefficient used to describe the height depen- 
dency of Tan in the surface mixing layer 


Standard deviation of the wind elevation angle at the nth 


point on the trajectory 


Power-law coefficient used to describe the height depend- 
ing of op, , in the surface mixing layer 
’ 


Power-law coefficient used to describe the height depen- 
dency of wind speed in the surface mixing layer 


Fraction of material reflected at the surface 


Washout coefficient or fraction of material removed by 
precipitation scavenging per unit time 


i ae ee te ae 
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Parameter 


TABLE 3-2 (Continued) 


Definition 
(b) Source Inputs 
Coordinates of the source 


Source Strength 


Stack exit temperature (°K) for the continuous source 


Inner radius of the stack at the exit 


Stack exit velocity 


Settling velocity for material in a given size category 


Decay coefficient or fraction of material lost per unit 
time 
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“ey features of the logic sequence shown in Figure 3-2 are that: 


e Solar altitudes for the release point and time of request 
are calculated using an algorithm suggested by Woolf 


(1968) | 


® The program then seeks the most recent cloud and ceiling 
data from the nearest reporting station to the release 


point 


e If there are no stations reporting cloud and ceiling data 
within a 50-kilometer radius from the release point or 
the data is more than three-hours old, the program 


selects a default net radiation index from Table 3-3 


e When cloud data is available, the program calculates 
an initial net radiation index based on solar altitude 
(also see Table 3-3) and adjusts the index based on the 


time-of-day total cloud cover, and ceiling height 


e The net radiation index and the wind speeds from the 


mesoscale wind-field prediction model for the release 


point and points along the cloud trajectory are used to 
select turbulent input parameters and wind power-law 


coefficients for the diffusion models 


8 Values of the surface mixing depth Hin, n for each nth 


point along the cloud trajectory are obtained from the 


mesoscale wind~field solution using the interpolation 


vrocedure given by Equation (2-5) 
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Figure 3-2, 


Schematic diagram illustrating procedures for obtaining 
inputs to the diffusion models, 
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TABLE 3-3 
NET RADIATION INDICES 


Daytime 
Index 
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Table 3-4(a) shows values of the ten-minute standard deviations of the 
wind azimuth angle in degrees (c, {rt = 10 min}) and associated values of the 
power-law coefficient m for describing o. variation with height above the sur- 
face based on the study of Swanson and Cramer. Values of the standard deviation 


of the wind elevation angle in degrees and the power-law coefficient n describing | 


the variation of o, with height above the surface are shown in Figure 3-4(b). 


| 
| 
These values of o. were assigned under the assumption that (after Cramer, et al., 
1972) turbulence is isotropic for short averaging times and that the one-fifth power ] 
law can be used as a scale factor to adjust an for averaging times less than 10 
minutes. Thus, | 
| 
1/5 
2.5 

= o. {r= 10 in} (252) se 
On oO, { min} (55 (3-6) 
om, {t = 10 min 1/2. 99 (3-7) 
The values of the power-law coefficient n were assigned under the assumptions 
that: : 

e Under stable conditions (NR = -1, -2) 


0 ~o Z 
ye E, ZR ts) (3-8) 

Z 

R 
j e Under near-neutral conditions (NR = 0,1) 
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TABLE 3-4 


TURBULENCE PARAMETERS FOR 
TRANSPORT/DIF FUSION MODELS 


(a) Ten-Minute Standard Deviations of the Wind Azimuth Angle in Degrees 
(a4) and Corresponding Values of the Power-Law Exponent m 


Net Radiation Index 


Wind Speed Be 
at 5 Meters 
(m/sec) 


20° =~. 109 A Bis, 9.8.6 =. 189 
20 =. L092 2s } » 12% ~12018.6 =.189)8.6 -.189 


LO Sade2 


14 -.105 


= LOO 


7 -.126|6.7 -. 287 
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(b) Standard Deviation of the Wind Elevation Angle In Degrees (o,,) and 
Corresponding Values of the Power-Law Exponent n i 


Wind Speed Net Radiation Index 
at 5 Meters 


(m/sec) 


EINER Th ce 


~ 


C 


(3-9) 


y o 
BE, Zz E, ZR 


e Under unstable conditions (NR = 2, 3, 4) 


Zz (3~10) 


The values of the wind power-law coefficient ae as a function of net radia- 
Y 
tion index and mean layer wind speed used in the program are shown in Table 3-5. 
These values are similar to those suggested for use at WSMR by Swanson and Hoidale 


(1962) and for Dugway Proving Ground by Cramer, et al., (1972). 


The plume-rise models described in Section 2,2.2 above require estimates 
of the ambient air temperature and, under stable conditions, the lapse-rate of poten- 
tial temperature. If surface meteorological data are available, the program auto- 
matically seiects the ambient air temperatures according to the logic shown in 
Figure 3-2 and the lapse-rate of potential temperature is calculated from the rawin- 
sonde sounding established using the logic described in Figure 3-1. If no surface 
data is available, the default values of surface temperature shown in Table 3~6 are 
used. The temperatures in the table are based on the monthly normal and maxi- 
mum temperatures at El Paso extracted from the report "Climatic Atlas of the 
United States'' published by the U. S. Department of Commerce in June 1968. The 
values shown for a net radiation index of -1 and -2 correspond to the normal mini- 
mum temperature for each month, Similarly, the values shown for net radiation | 
indices of 0 and 1 and for indices of 2 and 3 respectively correspond to the normal | 
average and normal maximum temperatures for each month, Finally, the surface 


temperatures for a net radiation index of -4 were obtained by taking an average of 
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TABLE 3-5 


WIND POWER-LAW EXPONENT p FOR 
TRANSPORT/DIF FUSION MODELS 


Mean Layer 
Wind Speed 
u (m/sec) 


Net Radiation Index 
aa ae ee a a ae 
<a sie 20 . 20 . 20 . 20 ad 


. 20 - 20 25 
kl! 
~L5 


a 


an eS ee 


Sapte mee eh feels Tr eae 


the recorded maximum temperature and the normal maximum temperature at El 


Paso for each month. 


If no rawinsonde data are available, the program selects default values of 
the lapse rate of potential temperature from Table 3-7 based on the net radiation 
index and the calculated wind speed at stack height. The potential temperature 
lapse rates shown in Table 3-6 are based on our experience in applying the Briggs 
stable plume rise formula in predicting plume rise from industrial stacks located 
in relatively warm and dry climates. No values of the default potential temperature 
lapse rates are given for wind speeds greater than 5 meters per second and for net 


radiation indices greater than 1 because plume rise is calculated from Equation 


(2-29) rather than Equation (2-32) under these conditions. 
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TABLE 3-6 
DEFAULT VALUES OF SURFACE TEMPERATURE Tg (OK) 


| Net Radiation Index 
| 


Month 
1, 0 | mh, a2 
Jan 292 286 279 272 
Feb | 296 290 282 275 
Mar | 300 294 | 286 277 
Apr | 303 299 | 290 282 
May | 308 | 304 | 295 287 
Jun | 312 308 300 292 
Jul | 312 308 301 294 
Aug 310 | 307 300 293 
Sep | 308 | 304 | 297 289 
Oct 303 299 291 283 
Nov | 297 292 | 284 | 275 
Dec 293 287 280 272 
| 
] 
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TABLE 3-7 


DEFAULT VALUES OF THE LAPSE RATE OF 
POTENTIAL TEMPERATURE # (deg K/m) 


Wind Speed Net Radiation Index 
ah} 


tm see} ! 0 -1 -2 
0~1.4 0.005 0,015 0, 030 0.040 
1.5 -2.9 0.003 0.010 0. 020 0.030 


3.0-4.9 0. 005 0.015 0.020 
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SECTION 4 


ANALYSIS, APPLICATION AND GRAPHICS 
PHASES OF EPAMS 


Figure 4-1 shows the executive components, in the form of a block diagram, 
of the EPAMS System. The EPAMS executive system manager, labeled XECAMS 4 
in Figure 4-1, is where all processing begins and ends. The executive systems man- 


ager on option passes control to the supervisors labeled ANALIZ, APPLY and 


PLTPRG which were designed and developed by the H. E, Cramer Company, Inc. 


to accomplish the technica! objectives of the contract. 


The analysis supervisor ANALIZ controls all processing in the Analysis 
Phase of EPAMS and passes control to the mixing layer and wind-field analysis 
foreman MIXLYR or to other meteorological analysis routines. The operational 


elements of MIXLYR are described in Section 4.1 below. The application super- 


visor APPLY passes control on option to the transport/diffusion application foreman 


TRYDIF or to other application routines. The operational elements of TRYDIF are 


described in Section 4,2. At the completion of either the Analysis or Application 
Phases, or at the completion of both phases, program control is returned to XECAMS., 
The executive systems manager can then address, on option, the graphics super- 


visor PLTPRG to obtain either printed or plotted solutions, or both types of solu- 


tions, developed in the application routines. Further details of the graphics super- 


visor PLTPRG are described in Section 4.3 below. 


4.1 ANALYSIS PHASE OF EPAMS 


As shown in Figure 4-1, the foreman MIXLYR of the Analysis Phase of 


7 


EPAMS is comprised of the two operational elements MXOBAN and SOLVMX. The 
operational subroutines of MXOBAN are shown in the form of a block diagram in 


Figure 4-2, The operational element MXOBAN controls the analysis of all meteor- 


Cee ee et 
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ological data and passes the results of the data analysis to the subroutine SOLVMX 
at the completion of execution. The EPAMS System user requests transport and 
diffusion information for a given time and source, which is then passed to MXOBAN 
via XECAMS, The MXOBAN mixing-layer analysis routine then proceeds to analyze 
rawinsonde, GWC and surface meteorological data according to the logic process 
described in Section 3.1. The subroutine MXUADT obtains rawinsonde data from 
the EPAMS data base for the desired solution time and source location. GWC pre- 
diction data and surface weather observations for use in the analysis are obtained 
from the data base via the subroutines MXGWDT and MXSFDT. These subroutines 
also use the GWC and surface data to modify the rawinsonde data if appropriate. 

The subroutine MXSFDT also calculates the net-radiation index for the surface met- 
eorological station using information obtained from the LONLAT and DAYNIT sub- 
routines. The rawinsonde data is then used in the subruutine LYRAHM to calculate 
the depth of the surface mixing layer and the mean wind components in the layer as 
described in Section 3.1. If insufficient data are available for calculating the 
mixing-layer depth or wind components, default data is supplied via the LYRDHM and 
MXSFDT subroutines. Control is then returned to the operational element MIXLYR 


and processing continues by branching to the subroutine SOLVMX, 


A block diagram showing operational subelements of SOLVMX is given in 
Figure 4-3, Processing begins through examination of the mesoscale model] input 
parameters passed to SOLVMX from MIXLYR to determine which of two procedures 
are to be used in initializing the mesoscale model, Thus, the input parameters from 
MIXLYR are compared in the subroutine MXSTRT with the input parameters used to 
generate the last solution of the mesoscale model resident in the EPAMS data base. 
The mesoscale model is initialized with the new inputs from MXOBAN if the mean 
layer wind speeds and elevation of the mixing depth above mean sea level (MSL) 
differ by more than twenty-five percent and the mean wind directions differ by more 
than ten degrees (cold start). The initial layer depth and momentum components for 


a cold start are given by the expressions 
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H = H' -h! (4-1) 


m3i, j n i,j 
= Ui he 4-2 
ae ( n ioe) ( ) 
A ee Sa ) (4-3) 
Ls n\n min 
where 
Ht = mixing depth (MSL) calculated for the present hour 


= terrain elevation (MSL) at each solution matrix grid 


ij 
“J point 
u = easterly component of the mean layer wind for the 
present hour 
h! _ — minimum terrain elevation (MSL) in the solution 
matrix 
We = northerly component of the mean layer wind for the 


present hour 


! If the differences are less than those specified above, subroutine MXSTRT calculates 


initialization parameters for a "hot" start from the expressions 


= ' we byt fs Aug Pm } t pa 
H (Hi hi) Bim He (4-4) 


(4-5) 


(4-6) 


whe re 
He . = mixing depth (MSL) at each solution matrix grid point 
xi for the previous wind-field solution 
H" = mixing depth (MSL) used to initialize the previous 
wind-field solution 
U. . — easterly component of the mean layer wind at each 
Vs solution matrix grid point for the previous wind- 
field solution 
U. = easterly component of the mean layer wind used to 
initialize the previous wind-field solution 
Vi . ~ northerly component of the mean layer wind at each 
a“ solution matrix grid point for the previous wind- 
field solution 
Vv = northerly component of the mean layer wind used 


to initialize the previous wird-field solution 


= = =" == 
7 


After the initialization procedure has been completed, the input parameters are 
passed to the subroutine MXSOLV via SOLVMX for execution of the mesoscale 


wind-field model described in Section 2.1. 
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subroutine MXSOLV will allow the mesoscale model to execute for 
a maximum model time of approximately 7200 seconds for a cold start and 
1800 seconds for a hot start or until the solution converges, whichever occurs 


first. The subroutine begins to check convergence criteria after 1200 seconds 


model time has elapsed for a cold start or 600 seconds for a hot start. The con- 
vergence criteria is then checked after every application of the nine-point filter 
(after every seventh time-step if the default value is used), Subroutine MXSOLV 


begins the convergence check by calculating the relative variance of the change in 


| 
mixing layer depth VR (AH, ) over the solution matrix according to the expres- | 
sion : 
\ nm ™! : 2 
V I I Elite eS 
R a an ) 2 Zs (Hast, j a 
IAL i 
ra $s 12) 
-nem BH ..=8 — : AS 
Pe 2 ( mi, j Heat \ (4-7) 


(n m n m 2| 
S(t 2 m| 3 Be J 
ter jar SST? SS ee 
j 
! 
where 
Hi = layer depth calculated for the present time step at each 
point in the n by m matrix 
j Hn = layer depth for the previous time step at cach point in the 


n by m matrix 


The subroutine determines the slope of the curve describing Ve (AH) as a func- 


tion of time and the difference between VR (HL) at the present and previous appli- 


Pe Rt RR ete 


cation of the filter. The solution is assumed to have converged if the 


j 
61 


ao UR ee 8 ee 


BS PR Ae 


onan 


lt ODL A A he 


ites ttateet salaries eee nineties een neeeiinimtimenaminaetaian Ailsa tb Nan ay yam. 


) Slope of the curve describing V,(44,,) as a function of 


time is equal to or less than zero 


e Present value of VR (au ) is less than the previous value 
m 
and the difference between the values is less than or equal 


to 1x 10° 


The program saves the converged mesoscale wind field solution and other param- 
eters calculated in the mixing-layer routines by passing them to the EPAMS data 
base. At this point the Analysis Routine is terminated and control is returned to 


XECAMS, 
4.2 APPLICATION PHASE OF EPAMS 


Inspection of Figure 4-1 above shows the foreman TRYDIF of the Applica- 
tions Phase of EPAMS is comprised of the two operational elements TRNPRT and 
DIFFUZ. These operational elements perform the calculations outlined in the block 
diagram shown in Figure 4-4 after receiving program control through the foreman 


TRYDIF from XECAMS and the supervisor APPLY. 


The subroutine TRNPRT calculates the trajectory of the cloud centroid 
downwind from the source location using the output from the mesoscale wind-field 
model and the procedures outlined in Section 2,2,1 above, The program calculates 
a point on the trajectory at intervals equal to one-tenth the spacing of the solution 
matrix of the mesoscale model (500 meters in the present case). At each of the 
calculated points on the trajectory, the program also determines the appropriate 
mixing-layer depth and wind speed for use in the diffusion models using the interpo- 
lation techniques described in Section 2.2.1. After completion of the calculations, 


program control is passed to subroutine DIFFUZ via TRYDIF. 
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The subroutine DIF FUZ calculates concentration and dosage with optional 
cloud depletion due to decay and either deposition losses from gravitational settling 
or precipitation scavenging using the instantaneous and continuous source models 
described in Section 2.2.2 above, At the users option, the program also calculates 
concentration and/or dosage isopleth half-widths at each point along the trajectory 
for instantaneous sources and concentration isopleth half-widths for continuous 
sources, At the completion of the diffusion model calculations the program writes 


the solutions to the EPAMS print file and returns control to XECAMS, 


4.3 GRAPHICS PHASE OF EPAMS 


Figure 4-5 shows a block diagram of the Graphics Phase of EPAMS. The 
program routine PLTPRG is the supervisor that controls all graphics output. The 
program retrieves solutions created in the application and Analysis Phases of 
EPAMS and prints and/or plots the solutions on option as directed by the user 
through XECAMS,. Options include the printing and plotting of intermediate and 
final wind-field patterns produced by the mesoscale model, isopleths of concentra- 
tion and dosage, cloud trajectories and mixing-layer depth contours. The plotted 
solutions are shown on a base map of the solution matrix grid with a background of 
terrain~height contours. The subroutine labled HEC004 in Figure 4-5 plots and 
labels the axes of each plot and the terrain height contours are added to the plot 
by subroutines HEC005 and HEC009, Wind velocity and direction vectors are cal- 
culated and plotted by subroutine HEC006. The program can, in subroutine HEC008, 


recalculate (at the user's option) additional cloud trajectories for other source loca- 


tions using the solution wind-field in the EPAMS data base. The supervisor PLTPRG 


can also calculate and display concentration and dosage isopleths for additional levels 
at the user's option. This feature is useful if the initially requested levels do not 
show details in which the user is interested or if the initial levels were either so 


large or so small that they do not occur in the solution matrix. It should be noted 
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that PLTPRG or its subroutines do not recalculate the wind-field nor the diffusion 
parameters, but use information from the last solutions to recalculate additional 


trajectories and isopleths. 


The subroutine HEC007 controls the plotting of the relative variance of 
the change in layer depth, the relative variance of the change in momentum, the 
total mass and total momentum as a function of model solution time. As noted in 
Section 4.1, the relative variance of the change in layer depth is used as criteria 
for establishing that the mesoscale wind-field solution has converged. Since the 
mesoscale wind-field model theoretically conserves mass and momentum, the plots 
of total mass and momentum versus model time are useful in monitoring the effects 


of numerical approximations in calculations of the wind-field solutions, 


At the completion of the plotting and print operations, program control is 


returned from PLTPRG to XECAMS. 
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SECTION 5 


EXAMPLE CALCULATIONS 


The Analysis and Application Phases of EPAMS were used to produce con- 
centration and dosage isoplith plots for four example cases. The calculations were 
made for meteorological conditions that occurred in the vicinity of WSMR during 


the 14th and 15th of November 1974, 


The surface weather pattern in the western United States on the 12th of 


November was characterized by a high-pressure ridge extending from the Pacific 
Northwest southeastward to Texas. The high-pressure ridge began to break down 
on the 13th of November with the rapid approach of a low-pressure system moving 
into the central U.S. from western Canada. The cold front associated with this 
system tended to remain east of the Rocky Mountains as it pushed south and east 
during the day. By the morning of the 14th, the low-pressure cell had moved into 
the northeast and the cold front extended southwestern through the eastern U.S., 
into Texas and northeast Mexico, then abruptly northward along the eastern slopes 
of the Rockies. The cold front had actually passed El Paso during the night from 
the "back-door", or moving from east to west. The passage was not accompanied 
by significant weather at WSMR, but broken middle and high clouds were reported. 
The high-pressure system behind the cold front had moved over WSMR by the mor- 
| ning of the 15th of November ani] calm winds were reported at many surface stations 


in the area. 


The temperature profiles ‘rom rawinsonde observations made at El Paso 
during the period are shown in Figure 4-1. Inspection of Figure 4-1 shown that 
the cold front had not reached El Paso by 1700 L on the 13th of November. By 

; 0500 Lon the 14th, the cold front had passed El Paso and the thermal inversion 
oecurring at about 800 millibars and associated with the frontal surface persisted | 


in the rawinsonde observations made at 1700 L on the 14th and 0500 L on the 15th : 


of November. Surface observations made at El Paso indicate that relatively shallow 


ote 


ee 


surface-based inversions formed during the nighttime and were eliminated by sur- 


face heating during the early morning hours, 


Three of the example cases were calculated for 0700 L, 0900 L and 1300 
L on the 14th of November to illustrate the performance of the Analysis Phase in 
analyzing and modifying the 0500 L rawinsonde observation made at El Paso on 
the 14th and the resulting effects on the concentration and dosage isoplith patterns 
produced by the Applications Phase of EPAMS. ‘The fourth case was calculated for 
0300 L, 15 November to illustrate program performance for very stable light-wind 
conditions at the surface. These calculations were all made for a hypothetical 
point source located about 31 kilometers southeast of WSMR using the default values 
for point source parameters given in Appendix A except the source strength Q was 
set to 100 units. The major terrain features of the solution grid are identified in 


Figure 4-2. The results of the calculations are described below. 
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5.8 INSTANTANEOUS POINT SOURCE CALCULATIONS FOR 0700L, 0900L 
AND 1300L, 14 NOVEMBER 1974 


Program requests were made for dosage and concentration calculations 
downwind from the hypothetical instantaneous point source southeast of WSMR for 


0700L, 0900L and 1300L, 14 November 1974, 


Because the request for a solution at 0700L was made for a time less than 
three hours after the 0500L rawinsonde observation at El Paso, the program used 
the temperature profile for 0500L shown in Figure 5-1 without modification of the 
surface, 850 millibar and 700 millibar temperatures and winds to calculate inputs 
to the mesoscale model. The rawinsonde analysis routine (MIXLYR) ignored the 
surface-based inversion indicated by the 0500L rawinsonde temperature profile 
because it was less than 100 meters in depth. Since the elevated inversion with a 
base at 832 millibars meets the criteria set for determining inversion bases 
described in Section 3.1, the analysis routine defines the mixing depth as 1697 
meters MSL or 491 meters above ground-level (AGL), The height-weighted mean 
wind components determined from the rawinsonde data below this level were u 
equal to -4.78 m ee and v equal 2.30 m ace or a wind direction of 116 


-1 
degrees at4.8 msec . 


The mesoscale wind-field model solution for 0700L is shown in Figure 5-3. 
Inspection of Figure 5-3 shows that southeasterly winds dominate the flow pattern 
in the scuthern portions of the Tularosa Valley and the valley Jornada Dei Muerto 
west of the Organ Mountains. Since the highest elevations of the San Andres, Organ 
and Sacramento Mountains extend above the depth of the surface mixing layer, they 
form an effective barrier to the southeasterly flow and divert the wind flow north- 
ward in the northern portion of the Tularosa Valley. The mesoscale wind-field 
model has set the wind components to zero and reduced the mixing layer to 30 


meters above terrain height in the higher elevations, which has resulted in down- 
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Figure 5-2. Map showing major topographic features of the mesoscale wind 
- field solution grid. 
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Mesoscale wind~field model solution for 0700L, 14 
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November 1974, 


slope flow as evidenced by the flow pattern just east of the Organ Mountains, A 
typical leeside depression of the mixing depth has formed to the west of the San 
Andres and Organ Mountains due to hydraulic jump effects and layer thickening 


has occurred east of the mountains. 


The results of the transport/diffusion model calculations for 0700L are 


shown in Figure 5-4. The cloud trajectory is diverted northward by the downslope 
flow occurring in the vicinity of the Organ Mountains and follows the south-south- 
westerly flow through the Tularosa Valley, Concentration isopleths, indicated by 
the heavier lines in Figure 5-4, are drawn for I x io: Lx io” and1x10 ~ 
source units per cubic meter. The dosage isopleths indicated by the thinner lines 
are drawn for 1 x 10, Lae 10° and 1x ie" source units per cubic meter. The 
isopleths were calculated using the diffusion model parametric inputs described in 
Section 3.2 associated with a net radiation index of -1 since the hour was nighttime, 


where nighttime is defined in the program as extending from one hour before sun- 


set to one hour after sunrise, and the cloud cover exceeded five-tenths. 


The inputs to the mesoscale wind-field mode} for the 0900L request for 
transport/diffusion calculations were still based on the 0500L rawinsonde measure- 
ment at El Paso, but the rawinsonde observation was modified because the observa- 
tion was more than three-hours old at 0900L. The modifications made by the mix- 
ing-layer analysis routine are shown in Figure 5-5, The routine has inse rted the 
0900L surface temperature and winds from El Paso and the GWC predicted 700 
millibar temperature and winds into the rawinsonde observation and has removed 
the 0500L temperature and winds observed at the surface and 812, 832, 877 and 
700 millibar levels. As a result, the analysis routine calculates a mixing depth 
of 2099 meters MSL (893 meters AGL) based on the increase in virtual temperature 
which occurs at 792 millibars. The mean wind components calculated for the mix- 


, : -1 -1 A 
ing layer were u equal -4.87 m sec and v equal 1.89 msec , ora wind 
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Figure S=4, Concentration isopleth (heavy line) and dosage isopleths fot 0700L, 14 
November 1974, Codes for isopleth levels are identified just below 
the computer plot. 
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El Paso rawinsonde observation at 0500L, 14 November 1974 (solid line) 
and modifications made by the mixing-layer analysis routine for the 0900L 
request for a transport/diffusion calculation (broken line). 
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The mesoseale wind-field model solution for 0900L is shown in Figure 
5-6. Since the wind components input to the model at 0900L are nearly identical 
to those used as inputs for 0700L, the solution shown in Figure 5-6 resembles the 
solution obtained for 0700L except in the vicinity of the Organ Mountains, The 
deeper mixing-layer depths obtained in the 0900L solution do not restrict the flow 
over the Organ Mountains nearly as much as in the 0700L solution and the down- 


slope flow on the east and south sides of the mountains does not occur, 


Figure 5-7 shows isopleths of concentration and dosage calculated for the 
point source assuming a release occurred at 0900L. The cloud trajectory is no 
longer affected by the downslope flow that occurred at 0700L and begins to travel 
westward, is diverted slightly southward around the Organ Mountains and then 
proceeds west-northwesterly in the Jornada Del Muerto valley flow patiern,. The 
concentration and dosage isopleths, drawn for the same levels used in Figure 5-4 for 
0700L, have become slightly more broad in the crosswind direction because the stan- 
dard deviations of the azimuth wind angle are significantly larger for a net radiation 
index of 2 used in the 0900L calculations than for the index of -1 used in the 0700L 
calculations, However, had the calculations been carried-out beyond the confines 
of the solution matrix shown in Figures 5-4 and 5-7 the isopleths for 0900L would 
have terminated or closed at shorter downwind distances because of the deeper 


{ mixing depths along the trajectory and the increased cloud widths. 


The latest rawinsonde data available for performing concentration and dos- 
age calculations for a release at 1300L was still the 0500L El Paso rawinsonde. 

Since the rawinsonde data was more than six hours old and the criteria for modify- 
: ing the rawinsonde and GWC data was not met, the mixingdayer analysis routine 
selected a default value of 1997 m MSL (92 AGL) from Table 3-1 for the mixing 


depth and the surface wind of 141 degrees at 3 meters per second from the surface 


‘ 
observations at El Paso for 1300L. The wind-flow pattern for 1300L produced by 
the mesoscale model, shown in Figure 5-8, indicates southeasterly flow over the 
76 
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Mesoscale wind~field model solution for 0900L, 14 November 1974, 
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Figure 5-7. Concentration isopleth (heavy line) and dosage isopleths for 0900L, 14 


November 1974, Codes for isopleth levels are identified just below 
the computer plot. 
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solution grid except immediately upwind of the San Andres and Organ Mountains 
and in the vicinity of the Sacramento Mountains, As indicated in Figure 5-9 
showing the concentration and dosage isopleths for a point source release at L300L, 
the cloud trajectory is directed northwestward and intersects the San Andres Moun- 
tains. The isopleths begin to expand rapidly at this point because the diffusion 
model compensates for attempted increases in concentration (in this case caused 
by rapidly decreasing mixing layer depth and insufficient increase in wind speed 

to compensate for reduced mixing depth) by increasing cloud dimensions in the 
plane of the horizon, The isopleths expand in a more regular fashion after the 
cloud has traversed the San Andres Mountains and reached the valley Jornado De 


Muerto. 


5.2 INSTANTANEOUS POINT SOURCE CALCULATIONS FOR 030uL, 15 NOV- 
EMBER 1974 


A program request was made for transport/diffusion calculations down- 


wind from the point-source southeast of WSMR for 0300L, 15 November 1974. 


“Since the latest rawinsonde observation available for analysis by the mixing-layer 


routine was more than six hours old (1700L, 14 November) and the GWC prediction 

data was rejected on temperature criteria, the program selected a default value of 

30 meters for the surface mixing-layer depth from Table 3-1 and a default wind 

component from the 0300L surface wind reported at El Paso. The surface wind 

speed at El Paso was reported as calm at 0300L and the mesoscale wind field model 

was initialized with zero momentum. The flow-field produced by the mesoscale | 
model is shown in Fimure 5-10. As indicated in the figure, downslope flow is occurr- | 
ing everywhere, but is greater west of the Sacramento Mountains and in the region 
just east of San Augustine Pass separating the Organ and San Andres Mountains. The 
cloud, as indicated in Figure 5-11 by the isopleth patterns, begins to move south~ 


west vard for a few kilometers, then moves toward the low elevations in the Tularosa 
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Figure 5-9, Concentration isopleth (heavy line) and dosage isopleths for 1300L, 14 


November 1974, Codes for isopleth levels are identified just below 
the computer plot. 
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Figure 5-10, Mesoscale wind~field model solution for 0300L, 15 November 1974. 
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Concentration isopleth (heavy line) and dosage isopleths for 0300L, 15 
November 1974. Codes for isopleth levels are identified just below 
the computer plot, 
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Valley at the southern boundary of the solution matrix grid, The highest wind 
speed along the cloud trajectory of 0.86 m a and mixing depth of 79 meters 
occurs at the boundary. Because the release occurred at nighttime with more 
than five-tenths cloud cover, the diffusion parameters for a net radiation index of 
-1 were used in the transport/diffusion routine. The relatively small values of the 
diffusion parameters in combination with the restricted mixing depth and low wind 
speeds lead to the high concentrations indicated by the narrow isopleths in Figure 


5-11. 
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APPENDIX A 


USERS' INSTRUCTIONS FOR THE ANALYSIS, 
APPLICATION AND GRAPHICS 
PHASES OF EPAMS 


This appendix gives users' instructions for the three main phases of EPAMS, ,| 


The computer programs that comprise the many phases of EPAMS are written in 


FORTRAN V and are designed for use on a Univac 1108 computer, A complete list~ , 
ing of the EPAMS computer program is given in Appendix B. Also, a block diagram | 
of the main components of each phase of EPAMS, as well as, temporary executive 
routines is given in Section 4 of the main body of this report. The program has been | 
delivered to WSMR with three dummy executive routines, these routines are shown 
as XECAMS, ANALIZ and APPLY in Figure 4-1, Section 4. The final design of these 
routines is not the responsibility of the H. E. Cramer Company and as such they are 
only presented here for completness and to provide temporary drivers for the Analy- 


sis, Application and Graphics Phases, See Section A.7 below for data file numbers 


assigned by the temporary executive XECAMS. 


The input variables required in the temporary executive XECAMS are des~ 


cribed in Section A.1 below. Also, optional input variables for the Analysis, Applica- 
tion and Graphics Phases are given in Sections A. 2, A.3, and A.4, respectively. An 
| example execution card deck is given in Section A.5 and the input data method used 


in the Analysis, Application and Graphics Phases is given in Section A.6 


A.1 EPAMS EXECUTIVE INPUTS 


The temporary executive XECAMS requires data input via the EPAMS input 
message file to direct the executive to execute one of the three phases of EPAMS. 
The set of variables given below are read by XECAMS prior to the execution of any 


of the EPAMS phases. One data card is read for each phase and only those variables 
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required by the particular phase need be input. A summary table of the input vari- 


ables and their respective units and default values is given in Table A-1, 


FORTRAN CARD 
VARIABLE FORMAT COLUMNS ro DEFINITION 


TYPE 1 13 1-3 This integer variable is 
used to direct the execu- 
tive to execute one of the 
EPAMS Phases, If set 
equal to 2, the Analysis 
Phase is executed and 
XECAMS calls the Analy- 
sis Supervisor ANALIZ. 
If set equal to 1, the Ap- 
plication Phase is exe- 
cuted and XECAMS calls 
the Application Supervi- 
sor APPLY. If set equal 
to 4 the Graphics Phase is 
executed and XECAMS 
calls the Graphics Super- 
visor PLTPRG. (This 
variable is referenced 
only in the executive 
XECAMS. ) 


TYPE 2 13 4-6 This integer variable is 

{ used to execute alternate 

sections of the Analysis 

i and Application Phases. 

If set equal to 1 or 2 and 
TYPE1 equals 2, the mix- 
ing layer and mesoscale 

wind field programs of 

4 the Analysis Phase are 

executed by a call to the 

mixing layer and meso- 
scale wind field programs. 

If set equal to 1 and TYPE 

1 equals 1, the transport/ 

diffusion programs of the 
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TABLE A-1 


SUMMARY TABLE OF INPUT VARIABLES 
FOR THE EPAMS TEMPORARY 
EXECUTIVE XECAMS 


FORTRAN 
VARIABLE 


CARD 
COLUMNS 


DEFAULT 
VALUE 


er 
month 

MGNTH Month corse 13-15 >1 and <12 

tHe a 

peas ei" 


FORMAT LIMITS 


TYPE 3 


2 3 
43 1G) 
ao] 


JULIAN Julian Day 


acti ss me — ioe 
HGUR Hour (GMT) = 200 and <23 N/A 


XPLACE 


2300000 and 
<450000 


Meters 360000 


23520000 and 
<3670000 


YPLACE 


ZPLACE 


Meters 


3585000 


Meters 


A-3 


‘ AD“-A038 399 CRAMER (H E) CO INC SALT LAKE CITY UTAH F/6 4/2 
MIXING“LAYER ANALYSIS ROUTINE AND TRANSPORT/DIFFUSION APPLICATI=-=-ETC(U) 
MAR 77 R K DUMBAULDr J R BJORKLUND DAAD07-76-C=-0023 
UNCLASSIFIED TR-76-106-01 ECOM-77=2 NL 
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FORTRAN CARD 
VARIABLE FORMAT COLUMNS DEFINITION 


TYPE2 Application Phase are exe- 
(Continued) cuted by a call to program 
TRYDIF, TYPE2 is not 
referenced in the Graphics 
Phase, 


| 
TYPES 13 7~9 This integer variable is | 
used to control input data | 
to parts of the Analysis, 
Application and Graphics 
Phases. If TYPE3 is set | 
less than or equal to zero, 
no input data is read by the | 
particular phase being exe- 
cuted, The programs will 
default or calculate all 
potential input data vari- 
ables. 


Specifically if TYPE3 is 
positive (greater than 

zero) and the mixing layer/ 
mesoscale section of the 
Analysis Phase is being 
executed, the data vari- 
ables given in Section A. 2 
are input via EPAMS file 
XINFIL (data card input). 
If the absolute value of 
TYPE3 is equal to 9, the 
program only executes the 
mixing layer solution. This 
option is used only to dis- 
play meteorological data 
and no solution files are 
generated or altered. 


If TYPE3 is positive 
(greater than zero), and 
the transport/diffusion 
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FORTRAN CARD 


TYPE3 
(Continued) 


Cen nn a ee ee ee 


VARIABLE FORMAT COLUMNS 


DEFINITION 


section of the Application 
Phase is being executed, 
the data variables given 
in Section A.3 are input 
via EPAMS file XINFIL 
(data card). If TYPE3 

is zero, the program 
defaults to an instantan- 
eous source with no input 
data. If equal to -1 the 
program defaults to a 
continuous SO» source 
given by the ASARCO 
copper stack in El Paso, 
Texas. If equal to -2, the 
program defaults to a con- 
tinuous SOg source given 
by the ASARCO lead stack 
in El Paso, Texas, 


If TYPE3 is positive 
(greater than zero) and the 
Graphics Phase is being 
executed, the variables 
given in Section A.4 are 
input via EPAMS file XIN- 
FIL (data card). 


If any of the three phases 
of EPAMS cannot complete 
its task, the particular 
phase will set TYPE3 equal 
to a negative number and 
return control to the execu- 
tive routines. If the phase 
solution was completed 
TYPES is greater than or 
equal to zero upon return. 
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FORTRAN CARD 
VARIABLE FORMAT COLUMNS DEFINITION 


JULIAN 13 10-12 The Julian day of the last | 
day of the month prior to 
the month of the desired | 
solution date and time. 
(Referenced only in Analy- 
sis Phase. ) 


a MGNTH 13 13-15 Month of the desired sol- 
ution. (Referenced only 
in Analysis Phase. ) 


DAY 13 16-18 Day of the desired solu- 
tion. (Referenced only 
i in Analysis Phase. ) 


H@UR 13 19-21 Hour of the desired solu- 
tion (past, present or | 
future) in Greenich mean 
time (00 to 23). (Refer- 
enced only in Analysis 
Phase. ) 


MINUTE 13 22-24 Minute of the desired sol- 
ution (00 to 60). (Refer- 
enced only in the Analysis 
Phase. ) 


XPLACE F10, 2 25-34 The UTM (Universal Trans- 
vers Mercator) easterly (x) 
coordinate of the source lo- 
cation in meters. If equal 
to blank or zero, the pro- | 
gram uses 360000 meters. 


a 


YPLACE F10.2 35-44 The UTM northerly (y) 
coordinate of the source 
location in meters. If 
equal to blank or zero, the 
program uses 3585000 
meters, 
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FORTRAN CARD 
VARIABLE FORMAT COLUMNS DEFINITION 


ZPLACE F10, 2 45-54 The elevation (MSL) of 
the source location in 
meters. If equal to 
blank or zero, the pro- | 
gram uses 1310 meters. 


A, 2 ANALYSIS PHASE INPUTS 


This section describes the variables input to the mixing layer part of the 
Analysis Phase, The variables given below are read only if the variable TYPE3 
from the executive is positive (greater than zero). These variables are read in the 
FORTRAN namelist format through namelist QLST1 via EPAMS file XINFIL (data 
card). The namelist format is explained in Section A.6 below. A summary table 


of the input variables is given with their respective units and default values in Table 


A-2, Possible input variables are: 


FORTRAN 
VARIABLE DEFINITION 
UN Vector component u of the mean layer wind for 


use in the mesoscale wind field calculation in 
units of meters/second. The component u is 
positive toward the east and negative toward the 
west. (Default = calculated from the meteorolo- 
gical data base.) 


UN = V-~> sin (02 180.0) 


whe re V = mean layer wind speed 
0 


= mean layer wind direction 
measured from 0 degrees 
north clockwise (direction 
from which wind is blow- 
ing) 
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TABLE A-2 


SUMMARY TABLE OF THE INPUT VARIABLES 
FOR THE ANALYSIS PHASE INPUT 
THROUGH NAMELIST QLST1 


INITIAL 
VALUE 


DEFAULT 
VALU 


Meters/ 

Second 

Mets rs/ 
So - 


FORTRAN 
VARIABLE 


UNITS LIMITS 


a ee 

= 
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(Table A-2 continued on following page. ) 


Q; the variable is equal to or less than the initial value the variable is set equal to 
the default value. 


Qrnese variables are normally calculated from the EPAMS data base, They are 
given as input variables to provide the capability to override calculated quantit- 
ies. 


TABLE A-2 (Continued) | 
FORTRAN 


VARIABLE LIMITS 


See previous page for footnote. 


DEFAULT 


id hk a a atl Ee 


TYPE 3=9 then 1 
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FORTRAN 


VARIABLE 


VN 


PN 


GN 


ABLKN 


SURF TN 


POTMPN 


ZELEVN 


DTLMDN 


DEFINITION 


Vector component v of the mean layer wind for 
use in the mesoscale wind field calculations in 
meters/second. The component v is positive 
toward the north and negative toward the south, 
(Default = calculated from the meteorological 
data base. ) 


VN = V-cos (9+ 180.0) 


Elevation (MSL) of the surface mixing layer above 
the meteorological station (if available) or source 
location in meters. (Default = calculated from 
the meteorological data base. ) 


Reduced gravity parameter for use in the meso- 
scale wind field calculations in units of meters/ 
second’, See Page 6, Section 2 of the main body 
of this report. (Default = calculated from the 
meteorological data base. ) 


Minimum surface mixing layer depth constraint 
in meters. This parameter must be greater than 
zero. (Default = 30 meters. ) 


Surface temperature in degrees kelvin. (Default = 
determined from the meteorological data base. ) 


Array containing the potential temperature at a 
maximum of 25 levels in degrees kelvin. (Default 
= determined from the meteorological data base. ) 


Array containing the elevation (MSL) in meters of 
each of the values input into POBTMPN. (Default = 
determined from the meteorological data base. ) 


Mesoscale model solution stability parameter 


shown in Equation 2-4 in the main body of this 
report. (Default = 0.4.) 
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FORTRAN 


VARIABLE 


ISMTHN 


NCNTN 


NETRDN 


KFLG1 


KFLG2 


DEFINITION 


Number of mesoscale model time iterations 
between applications of a nine-point low-pass 
filter. (Default = 7, filtered every 7th step.) 


Number of mesoscale model time iterations 
between recalculations of the model time incre- 
ment constant. (Default = 2, recalculated every 
2nd step.) 


Net radiation index as defined in Section 3, 2 in 
the main body of this report. (Default = calcu- 
lated from the meteorological data base, ) 


Switch to control the initialization of the meso- 
scale model. If set equal to zero, the program 
will initialize the mesoscale model with a cold 

or hot start according to the rules given in Sec- 
tion 4.1 of the main body of this report. If KFLG1 
is set equal to 1, the program forces a cold start, 
i.e., the mode} is initialized with the constant 
momentum field calculated from UN, VN and PN 
above. (Default = 0.) 


Switch to control the print output from the Analysis 
Phase. If set equal to 0, only completion mes- 
sages are output to the EPAMS output file XOTFIL 
(printer), If set non-zero, the Analysis Phase 
prints detailed tables of the meteorological data 
used and solution characteristics of the mesoscale 
model. If the absolute value of TYPE3 defined 

in Section A.1 is equal to 9, KF LG2 is automat- 
ically set non-zero. However, no mesoscale out- 
put is produced. (Default = 0.) 


APPLICATION PHASE INPUTS 


This section describes the variables that are input to the transport/diffusion 


section of the Application Phase. The variables given below are read only if TYPE3 
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from the executive is positive (greater than zero). These variables are read in 


the FORTRAN namelist format through namelist QLST2 via EPAMS input file 


XINFIL. A summary table of the input variables is given with their respective 


default values in Table A-3, Possible input variables are: 


FORTRAN 


VARIABLE 


XS@ZUR 


YS@GUR 


ZSGUR 


SIGXG 


SIGYG 


DEFINITION 


UTM easterly (x) coordinate of the source location 

in meters. (Default = value used in Analysis Phase 
or if continuous source 1 or 2 is specified = 356000 
meters (ASARCO), ) 


UTM northerly (y) coordinate of the source location 
in meters. (Default = value used in Analysis Phase 
or if continuous source 1 or 2 is specified = 3520100 
meters (ASARCO), ) 


Elevation (MSL) of the surface at the source location 
in meters. (Default = value used in Analysis Phase 
or if continuous source 1 or 2 is specified = 1150.6 
meters (ASARCO).) 


Source strength of source. Input in units of total 
material released if an instantaneous source or 
total material per second if a continuous source. 
(Default = 1 if an instantaneous source, = 2520 
g/second SO5 if continuous source 1 (ASARCO 
Copper), = 1764 g/second SOg if continuous 
source 2 (ASARCO lead). ) 


Standard deviation of the alongwind source dis- 
tribution after cloud stabilization for instantan- 
eous sources only. Input in units of meters. 
(Default = 0, ) 


Standard deviation of the crosswind source dis- 
tribution after cloud stabilization for instantaneous 
sources only. Input in units of meters. (Default = 
0.) 
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TABLE A-3 


SUMMARY TABLE OF THE INPUT VARIABLES 
FOR THE APPLICATION PHASE INPUT 
THROUGH NAMELIST QLST2 


FORTRAN DEFAULT 


ITS LIM 
VARIABLE ici YT VALUK2) 


>300000 and 


XSGUR Meters <450000 


XPLACE 


23520000 and 


Mere <3670000 


YPLACE 


=1 if KINDSG = 0 
=2520 g/sec. if KINDSG = 1 
=1764 g/sec. if KINDS@ = 2 


=0 if KINDS@ = 0 
H Meters >0 =252.4 if KINDS@ = 1 
=186.5 if KINDSG = 2 


(Table A-3 continued on following page. ) 


Orne source strength input unit is total material released for instantaneous sources 
(KINDSG@ = 0) or total material released per second for continuous sources (KINDS@ 
> 0). If the source strength is input as grams for instantaneous sources or grams 
per second for continuous sources, the output concentration units are grams per 
cubic meter and dosage units are gram seconds per cubic meter. Dosage is cal- 
culated only for instantaneous sources, 


ees 


Q; the value of the variable is less than or equal to zero, the variable is set equal 
to the default value. 
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TABLE A-3 (Continued) 
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FORTRAN ca : DEFAULT 
VARIABLE si de Bugle VALUE 


=0 if TYPE3 > 0 
=1 if TYPE3 = -1 
=2 if TYPE3 = -2 


=366,3 if KINDSG = 1 
=352.4 if KINDS@ = 2 


=4, 877 if KINDS@ = 1 
AD Mete >0 
Meters/ =2.1476 if KINDSG = 1 
saat an Second =2,7512 if KINDSG = 2 
| 
N 
vs leters/ 
H Second 


i Concentration 1x 1074, 1x 1078 
CONISS Output Units and 1 x 10712 


STTEMP > Ambient Air 


Dosage 1x 1074, 1 x 1076 


S gy 
D@SISS Output Units and 1 x 10714 


i 5 (Table A-3 continued on following page. ) 
} Qe previous page for footnote. 
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{ TABLE A-3 (Continued) 


FORTRAN DEFAULT 
VARIABLE VALUE 


a 
Paver [ [wes 
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Parve [meee [moi fe 


Qeo previous page for footnote. 
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FORTRAN 
VARIABLE 


SIGZG 


H 


TAU 


ALPHA 
BETA 


KINDS@ 


DECAY 


BLAMDA 


Sec ibid gh inert ae ye 


DEFINITION 


Standard deviation of the vertical source distribu- 
tion after cloud stabilization for instantaneous 
sources only. Input in units of meters. (Default 
— 0: ) 


Source emission release height above ground if a 
continuous source or effective height above ground 
after cloud stabilization for instantaneous sources, 
Input in meters. (Default = 0 if instantaneous, = 
252.4 meters if continuous source 1 (ASARCO 
Copper), = 186.5 meters if continuous source 2 
(ASARCO Lead), ) 


Height in meters above ground at which concen- 
tration and/or dosage is to be calculated. (Default 
=" ())2 ) 


Source emission time in seconds for instantaneous 
sources only. (Default = 2.5 seconds. ) 


Crosswind diffusion coefficient. (Default = 1.) 
Vertical diffusion coefficient. (Default = 1.) 


Source type switch, where if KINDS@ = 0, the pro- 
gram assumes an instantaneous source. If KINDS@ 
is greater than zero, the program assumes a con- 
tinuous source. If set equal to 1, the program 
assumes a continuous source 1 which is the copper 
smelter stack at ASARCO, El Paso, Texas. If set 
equal to 2, the program assumes a continuous 
source which is the lead smelter stack at ASARCO 
El Paso, Texas, (Default = 0.) 


Coefficient of time dependent decay input as a frac- 
tion. (Default = 0.) 


Coefficient of precipitation scavenging input as a 
fraction. (Default = 0.) 
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FORTRAN 


VARIABLE 


GAMMA 


STTEMP 


RADIUS 


VS 


CA@NISS 


D@SISS 


TEMPSF 


RTLAPS 


KFLGL 


DEFINITION 


Coefficient of vertical reflection input as the frac- 
tion of material reflected at the surface. (Default 
= Le) 


Stack effluent exit temperature in degrees Kelvin 
for continuous sources only, (Default = 366.3 °K 

if continuous source 1, = 352.4 if continuous source 
2, no default is provided for KINDSG@ > 2.) 


Internal stack radius in meters for continuous 
sources only. (Default = 4.877 meters for con- 
tinuous source 1, = 3.658 meters for continuous 
source 2. No default is provided for KINDS@ > 2.) 


Particle or droplet settling velocity in meters/ 
second for concentration and/or dosage with deple- 
tion due to gravitational deposition. (Default = 0.) 


Array of a maximum of 3 concentration values for 
which concentration isopleth half widths are to be 
calculated, (Default = 1x 1074, 1x 10-6 and 

1x 10712.) 


Array of a maximum of 3 dosage values for which 
dosage isopleth half widths are calculated for 
instantaneous sources only, (Default = 1 x 1974, 
1x 1078 and 1 x 10712, ) 


Surface temperature in degrees Kelvin for con- 
tinuous source only. (Default = calculated in 
Analysis Phase or from historical data. ) 


Lapse rate of potential temperature in degrees 
Kelvin/meter for continuous sources only. 
(Default = calculated from eievation and potential 
temperature provided by the Analysis Phase or 
from historical data. ) 


Mesoscale model! solution switch. If set equal to 0, 


the program will calculate the source cloud trajec- 
tory from the present or most current mesoscale 
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FORTRAN 
VARIABLE DEFINITION 
KFLG1 wind field solution. If non-zero, the program will 
(Continued) use the past mesoscale solution. (Default = 0.) 
KFLG2 Detail print output switch. If equal to zero, the 
program writes to the EPAMS output file the calcu- 
lated concentration, dosage, standard deviation of 
the lateral concentration distribution and the isopleth 
half widths at each trajectory point. If non-zero, 
the program also outputs the trajectory points with 
the wind speed, surface layer depth, time and total 
distance. (Default = 0.) 
A.4 GRAPHICS PHASE INPUTS 


This section describes the variables that are input to the Graphics Phase. 
The variables given below are read only if TYPE3 from the executive is positive 
(greater than zero), These variables are read using the FORTRAN namelist for- 
mat through namelist QLST3 via EPAMS file XINFIL. A summary table of the 
input variables is given with their respective default values in Table A-4. The 
plot routines assume all plot output is to be written to FORTRAN logical unit 1 
on magnetic tape. Also, the plot routines assume the plotter is a drum plotter 
with at least 24 inch paper and with three different pen turrets to hold optional plot 


colors. The best color combination is for pen 1 to be black, pen 2 to be red and 


pen 3 to be blue. 


FORTRAN 
VARIABLE DEFINITION 
i 
; ISW Array of print and/or plot options 
ISW(1) If equal to zero, mesoscale wind field solutions i 


are not printed. If equal to 1, mesoscale wind 
field solutions are printed in wind speed vector 
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TABLE A-4 


SUMMARY TABLE OF THE INPUT VARIABLES 
FOR THE GRAPHICS PHASE INPUT 
THROUGH NAMELIST QLST3 


ee me LIMITS DEFAULT VALU 


IS\W(8) 


ISW(9) 


ISW(10) 


>300000 and 


MV 2 
Teters <450000 


(Table A-4 continued on following page. ) 


Decker contour intervals are defaulted to 1219.2, 1524, 1828.8, 2133.6, 2438.4, 
2743, 2, 3048, 3352, 8 meters. 


Ohrixing layer depth contour intervals are defaulted to 50, 100, 200, 300, 400, 600, 
800, 1000, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2600 meters, 


De the value of the variable is less than or equal to zero the variable is set equal 
to its default value. 


A-19 


“5 “Aiea meta 


AB SE NBT TN AN AE I IR AN AOR ORR 


Lian Ba 8p Ae? <= 


TABLE A-4 (Continued) 


FORTRAN 
i ‘S iy! no 1 , 
ieee UNIT LIMITS DEFAULT VALUE 
>3520000 and 


NCGNTS Integer 
Meters 
Inches/ Meter 


CNISS COBRENEESEIOS 
Units 
Dosage 
Units 


Qrerrain contour intervals are defaulted to 1219.2, 1524, 1828, 8, 2133.6, 2438, 4, 
2743.2, 3048, 3352.8 meters. 


Qhrixing layer depth contour intervals are defaulted to 50, 100, 200, 300, 400, 600, 
800, 1000, 1200, 1400, 1600, 1800, 2000, 2200, 2400, 2600 meters. 


OF the value of the variable is less than or equal to zero the variable is set equal 
to its default value. 
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FORTRAN 
VARIABLE DEFINITION 


ISW(1) components for each elapsed 30 minutes of meso- 

(Continued) scale model time plus the final solution. If equal 
to 2, the wind speed and direction are printed 
rather than vector components. If negative, only 
the final wind field solution is printed. (Allow for 
31 pages of output for each mesoscale wind field 
solution printed. The maximum is 155 pages. ) 
(Default = 0.) 


ISW (2) If non-zero, isopleths of concentration and/or 
dosage are plotted on a terrain contour map, The 
data is taken from the transport/diffusion solution 
file. Also, isopleths can be recalculated and plot- 
ted by the use of CNISS and/or DSISS below. (Default 
= al) 


ISW(3) If non-zero, the relative variance of the change in 
momentum and the relative variance of the change 
in layer depth are plotted as a function of mesoscale 
model solution time. For an explanation of the pur- 
poses of this plot see Sections 4.2 and 4.3 of the 
main body of this report. (Default = 0.) 


ISW (4) Not used at present. 


ISW(5) If non-zero, the total momentum and total volume 
from the mesoscale solution are plotted as a func- 
tion of mesoscale model solution time. For an 
explanation of the purpose of this plot, see Section 
4.3 of the main body of this document. (Default 
= 0.) 


ISW (6) Not used at present, 


ISW(7) If non-zero, the program uses elevation (MSL) 
rather than depth in plotting contours of the mixing 
layer under ISW(10) below. (Default = 0.) 


ISW(8) If non-zero, a contour plot of the terrain elevation 
(MSL) is generated. (Default = 0.) 
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FORTRAN 


VARIABLE 


ISW(9) 


ISW (10) 


ISW(11) 


ISW(12) 
through 
ISW(19) 


XST 


YST 


DEFINITION 


If equal to 1, wind speed and direction vectors are 
plotted on a terrain contour map. If equal to 2, | 
the wind speed and direction vectors are plotted | 
without the terrain contours. If ISW(9) is positive, | 
each intermediate mesoscale solution is plotted at 

30 minute intervals including the final solution. A 

maximum of 5 plots could be generated under this 

option. However, if ISW(9) is negative, only the 

final solution is plotted. If this option is selected, 

specify 3 different pen colors for the three differ~ 

ent pen turrets (black, red and blue). (Default = 

0.) 


If equal to 1, contours of layer depth or elevation 
depending on ISW(7) are plotted on plots generated 
by ISW(9) above. If equal to 2, contours of layer 

depth or elevation are plotted on contours of ter- 

rain elevation for each of the intermediate meso- 
scale solutions. If equal to -2, only the final sol- 
ution is plotted. (Default = 0.) 


If non-zero, trajectories generated from optional 
start locations given by XST and YST below are 
calculated and plotted on plots generated under 
ISW(9) above. If ISW(11) is positive, trajectories 
are plotted only. If negative, trajectories are both 
printed and plotted. (Default = 0.) 


Not used at present. 


Array of a maximum of 5 UTM easterly (x) coordin- 
ates of the start location of optional trajectories in 
meters. 


Array of a maximum of 5 UTM northerly (y) coor- 
dinates of the start location of optional trajectories 
in meters. 
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FORTRAN 
VARIABLE DEFINITION 


NCONT Number of terrain elevation contours to generate. 
(Default = 8.) 


CONT Array of a maximum of 8 contour elevation levels 
in meters. (Default = 1219.2 to 3352. 8 in incre- 
ments of 304. 8 meters or 1000 feet. ) 


NCGNTS Number of mixing-layer depth (elevation) contours 
to generate. (Default = 16.) 


C@NTS Array of a maximum of 16 layer depth (elevation) 
contour levels in meters. (Default = 50, 100, 200, 
300, 400 to 2600 in increments of 200 meters. ) 


SCL Plot scale factor for ISW(2), ISW(8), ISW(9), ISW 
(10) and ISW(11) in inches per meter. (Default = 
1.0 x 10-4, ) 

CNISS Array of a maximum of 3 concentration isopleth 


values to use in place of those in the diffusion 
solution file under ISW(2) above. (No defaults 
provided. ) 


DSISS Array of a maximum of 3 dosage isopleth values 
to use in place of those in the diffusion solution 
file under ISW(2) above. (No defaults provided. ) 


| A.5 EXAMPLE EPAMS EXECUTION CARD DECK 


Figure A-5 gives an example job card deck that executes the EPAMS pro- 


~ 


gram. The first three cards shown in Figure A-5 initialize the job and define the 


EPAMS program file, The next two cards are used to assign a plot output tape to 
be generated by the Graphics Phase of EPAMS. The sixth card shown assigns all 
of the catalogued data and solution files required by EPAMS and is an element 


residing in the EPAMS program file. The seventh card executes the data case 
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cards, The first two data cards direct the Analysis Phase to analyse the meteoro~ 
logical data on November 13, 1974 at 1400Z. Also, the source location and eleva- 
tion is given on data card one to specify the release point and enable the program 

to select the closest meteorological data, Optional input was selected for the Analy- 
sis Phase in order to set KF LG2 to 1 to obtain detailed output from this phase. The 
next two data cards direct EPAMS to execute, with detailed output, the Application 
Phase for the concentration and dosage from an instaneous source, ‘The data cards 
give a source strength of 100 and allow the program to default all other possible 
inputs. The last two data cards execute the Graphics Phase and specify plots of 
concentration and dosage isopleths are to be produced, as well as plots of the final 
mesoscale wind field solution, relative variance of momentum and depth and total 


momentum and volume, 


The order of the data cards in the example gives an EPAMS solution from 
start to end, However, the user could execute any one or group of the EPAMS 
phases depending on the solution desired, An example would he to place a set of 
Application and Graphics data cards after those shown in Figure 5-1 to produce 
concentration and dosage for an alternate source location. The user in this case 
would duplicate the last 4 data cards and would then include the alternate source 


location in the $QLST2 card of the second set. 
A.6 DATA INPUT METHOD 


The three phases of EPAMS use the FORTRAN namelist data input method. 


The namelist input data must be in a specific form in order to be read using a 


NAMELIST list. The first character in each card to be read must be blank, The 
first card in the namelist list contains the namelist name preceded by the character 
$ (SQLST1), The last card in each namelist list contains SEND to terminate the 


list. The form of the remaining data items in the list may be: 
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a. Variable Name = Constant - The variable name may be a sub- 
scripted array name or a single variable name, Subscripts must 


be integer constants. Theconstant may be integer or real. 


b. Array Name = Set of Constants (Separated by Commas) - The 
array name is not subscripted. The set of constants consists of 
constants of the type real or integer. The number of constants must 
be less than or equal to the array size. Successive occurrences of 


the same constant can be represent in the form k * constant. 


The sequence of the input data parameters within the list is not significant. 

A more detailed explanation of the FORTRAN NAMELIST can be found in any Fortran 
Language Manual. The input parameters within the namelist that have default values 
are initialized to zero or other values in the Analysis Phase (Table A-2) prior to 
input of the first case. Parameters that are not used or have default constant or 
calculated values need not appear in the namelist list. When multiple cases are exe- 
cuted, all parameters retain their values from the previous case except in the 
Analysis Phase and are changed only by input. Also, the computer program assumes 


that core storage has been zeroed prior to execution. 
Beal EPAMS TEMPORARY EXECUTIVE AND DATA FILES 


The H. E. Cramer Company has provided a temporary executive routine 
XECAMS to drive the programs that comprise the Analysis, Application and Graphics 
Phases of EPAMS. Program XECAMS is the main executive driver of EPAMS. This 
program assigns logical file numbers to each of the 12 EPAMS da.a and solution 
files. The internal file names are held in common block FILESdefined on page B-8 
of Appendix B and are assumed defined prior to the execution of any EPAMS phase, 


The internal file name, logical number and contents are 
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FILE NAME 


XINFIL 


XOTFIL 


OBFILE 


UAFILE 


SFCOBF 


SBFILE 


GWCPFL 


GWCFIL 


WSTRRN 


ao 


21 


10 


11 


12 


LOGICAL UNIT NUMBER 
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CONTENTS 


Card reader, 


Printer, 


Upper-air data master file directory. 
This file contains the station number, 
station location and hours for which 
upper air data has been recorded. 


Upper-air data inventory file contain- 
ing the observations for the stations 
and hours identified in QGBFILE. 


Surface data master file directory 
containing the station number, station 
location and hours for which surface 
data has been recorded. 


Surface data inventory file containing 
the observation for the stations and 
hours identified in SFC@BF. 


GWC (Global Weather Central) data 
master file directory containing the 
location and rawinsonde (upper-air) 
hours for which GWC prediction data 
has been calculated. 


GWC data inventory file containing the 
prediction data for the 18 hours follow- 
ing each rawinsonde hour and for each 
location given in GWCPFL. 


White Sands Missile Range terrain 
elevation (MSL meters) data, This 
file contains the 41 x 41 (1681) point 
matrix of elevation data where the 
internal 31 x 31points are real terrain 
data and the 5 outermost rows and 
columns are artificially filled with the 
minimum WSMR elevation. The data 
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FILE NAME 


MXLYRF 


TRDFFL 


14 


LOGICAL UNIT NUMBER 


A-28 


CONTENTS 


are stored starting with the most 
south-west corner and proceeding 
eastward, then moving one point 

north and repeating so the last point 

is the most north-east elevation. 

The real terrain data of the internal 

31 x 31 points has a south-west UTM 
coordinate of 300000, 3520000 meters. 
The last of the real terrain to the 

east is at x = 450000 and to the north 
is y = 3670000 covering 150 kilometers 
in both directions. This file is treated 
as a read-only file by the EPAMS 
phases. 


This file contains the present and past 
solutions of the mesoscale wind-field 
model. The file contains two data 
blocks, both 5111 words long. Block 1 
is the past solution and block 2 is the 
present or most current solution. Each 
block consists of the variables given in 
common block MXBUFR defined in the 
FORTRAN procedure listing in Appendix 
B, page B-9, 


This file contains the diffusion model 
solutions from the Application Phase. 

The file contains one block of 4090 words. 
The block contains the variables NCALC, 
IULIAN, IMQNTH, IDAY, IHGUR, 
NINUTE, XSGURC, YSG@URC, ZS@URC, 
KINDS, CQNI (3), D@SI(3), DUMMY (664), 
XX (310), YY(310), SIGMY(310), CHIT 
(310), D@ST(310), YCHIT(310, 3) and 
YDUST(310,3). The variable NCALC is 
the total number of points in the cloud 
trajectory, TULIAN is the Julian day of 

the solution, IMONTH is the month, IDAY 
is the day, THOUR is the hour, NINUTE is 
the minute and XS@URC, YS@URC and 
ASQAUR is the source location. Also, 


be ecmteyed 


FILE NAME LOGICAL UNIT NUMBER CONTENTS 


KINDS is the source type number, 
CQYNI concentration isopleth levels, 
DQSI dosage isopleth levels and 
DUMMY is extra space for future use. 
XX and YY define the trajectory loca- 
tion, SIGMY is the standard deviation 
of the crosswind concentration (dosage) 
distribution at XX, YY, CHIT and 
DYST are the concentration and dosage 
and YCHIT and YDQST are the con- 
centration and dosage isopleth halfwidths, 
respectively. 


SPARE 9 This file is used to store intermediate 
output from the mesoscale wind-field 
mode! for subsequent printing and/or 
plotting by the Graphics Phase. The 
file can contain from 2 to 6 blocks of 
data. The last block is indicated by 
the first word in the block being negative. 
This block contains the relative variance 
of the change in momentum and depth for 
each model time step, as well as the 
total volume and momentum for each 
time step. Each block contains 5111 
words where the first blocks are the 
variables in common MXBUFR defined 
on page B-9 of Appendix B. The last 
block is also the same as common 
MXBUFR except the common block 
SWELL defined in subroutine MYPRG 
(page B88, Appendix B) has replaced the 
arrays M, N and P in MXBUFR and the 
first word of the block is set to the nega- 
tive of the total number of mesoscale 
model time steps. Also, the first blocks 
are written at approximately 30 minutes 
model time intervals until model conver- 
gence. The final solution is then written 
before the last data block of this file. 
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Included in the common block of internal file names are three variables that define 


the number of data blocks in the master directory files. These variables are 
MXBLKU, MXBLKG and MXBLKS for files GBFILE, GWCPFL and SFCQOBF, re- 
spectively. The Analysis Phase programs use these variables to start reading 

the master file directories with the most current data, where the most current data 
is assumed to be at the end of the files. In a real time environment, these varia- 
bles would probably be set to 1 with the most current observations being the only 
data in the files. These variables are also initialized in the executive EXAMS on 


page B-22 of Appendix B. 


A-30 


Re 


APPENDIX B 


Appendix B contains the complete FORTRAN listing of the EPAMS computer 
programs. Page B-2 shows the UNIVAC 1108 element (FILEPAMS) used to assign 
the necessary data files. Page B-3 and B-4 show the UNIVAC 1108 element (EPAMS- 
MAP/WSMR) used to collect the EPAMS programs into a segmented absolute pro- 
gram named ABSEPAMS,. Pages B-5 through B-21 show the UNIVAC 1108 FORTRAN 
procedure used to define the common blocks used by the EPAMS programs, Also, 
the procedure defines all important variables used in the program, The remainder 
of the Appendix contains a FORTRAN listing of each program used in the EPAMS 
system. Also, the subroutines UADATA, XTRCUA, SATVP, SFCDTA, TESTFD, 
GWCDTA, BLKIN, BLKQGUT and DRMRD are routines written and supplied by per- 
sonnel at WSMR. However, some of these routines have been modified and are not 


interchangeable with previous copies of the routines. 
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